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Abstract. We study the dynamical properties of small regulatory networks treated as 
non autonomous dynamical systems called modules when working inside larger networks 
or, equivalently when subject to external signal inputs. Particular emphasis is put on the 
interplay between the internal properties of the open systems and the different possible inputs 
on them to deduce new functionalities of the modules. We use discrete-time, piecewise-affine 
and piecewise-contracting models with interactions of a regulatory nature to perform our 
study. 
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1. Introduction 

The structure of genetic regulatory networks can be abstracted by directed graphs, where the 
nodes represent genes and the arrows (oriented edges) stand for their interactions through 
transcription/translation products. These interactions may be either activations or inhibi- 
tions. Given the large number of components in most networks of biological interest, con- 
nected by positive and negative feedback loops, the comprehension of the dynamics of a 
system is often difficult if not impossible. In this context, mathematical modeling eventually 
supported by computer tools can contribute to the analysis of a regulatory network by allow- 
ing the biologist to focus on a restricted number of plausible hypotheses, or to easily read 
some observed features of the system. We refer to |2| j8] and references therein for an account 
of the huge and still growing literature on this topic. 

It has been evidenced that the existence of network motifs, a set of recurring patterns in- 
side large biological regulatory networks, shall give new insight on the understanding of the 
performances of the global network. We refer to the review [3] (see also [2]) for a clear and 
complete description of such a point of view and references therein for recent work in differ- 
ent specific biological contexts. Since most of the time these modules [19] also have input 
and output interactions relating to the rest of the network, they may be considered as open 
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dynamical systems, meaning distinguished parts of larger 'closed' systems, the remaining of 
the network being an 'environment'. 

In the classical literature of Dynamical Systems, "open" or "forced" dynamical systems are 
defined as those for which the dynamics depends on internal rules as well as on inputs from the 
environement, then producing some outputs back to the environement. The notion of open 
systems is also the groundwork of System Theory and many results concerning the possibility 
to control linear and nonlinear systems have been obtained in the last decades [20} LT3] . Here 
we are interested in a description of the dynamics of these systems in all the parameter space 
and subjected to any external signal. 

When comparing these modules with the standard corresponding regulatory networks we are 
faced to new questions to understand their dynamical behavior and, therefore, their possible 
new functionalities. Here we study some of these questions by searching first for general 
properties of such systems and by trying to identify the fundamental mechanisms at work in 
the open systems we deal with. Although our definition of module is not the most general 
possible, as it will appear next, it seems that it contains all the cases treated in the literature 
and, remarkably, it allows simpler proofs in many cases. 

These general properties must be complemented in each particular case with further analysis 
to derive additional specific properties of a given module and we perform this task for some 
examples. 

What we found in the general case may be summarized as follows: each elementary input 
selects a particular subset of the phase space and a particular dynamical rule among a (finite) 
set of possibilities allowed by the module; therefore any (finite or infinite) sequence of inputs 
drives the dynamics of the module along a pre-determined maze of possible paths. Simple 
sequences of inputs give rise to simple dynamics of the module and therefore to simple 
responses (outputs), but, since they fix the system in a "corner" of its phase space, the 
system may acquire new functionalities as we shall see in some examples. Moreover, because 
such constraints in phase space may depend of the values of the parameters of the system 
and not only of its structure (the pattern of interactions) it turns out that the structure 
cannot determine the functionality in all cases [12]. It will also be clear from what follows 
that, on the contrary, in some cases it can happen that the type of response is not essentially 
affected by a change of parameters, and only particular aspects do, e. g. speeding or slowing 
of the response. Interestingly, when a module can perform more than one function, each one 
is generally robustly provided. This fact does not exclude the possibility of a fine tuning of 
the interaction between the input and the internal dynamics that may end in an interesting 
interplay between both. 

Despite the fact that the proofs of our results is done for a special type of model of regulatory 
networks, i. e. discrete time dynamical systems [HJ El [HI [25] , we believe that these results 
have a more general thrust, for the underlying mechanisms are common to other types of 
dynamical systems, as ODEs for instance. 

Besides, discrete-time models provide a simple framework where the consequences of inter- 
action delays are already included and they largely benefit from tools and techniques in the 
Dynamical Systems theory [8[ [6] . 
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We study two different types of modules. The first one, called regulatory cascades (RC), 
is defined as any open network without (internal) circuits. The second type, called forced 
circuits (FC), is any open circuit. The general case can be build as a combination of such 
elementary pieces and at the end we sketch the analysis of one such mixed RC-FC module 
on the bases of our previous results. 

For the RC modules we give a complete description of the dynamics. We specify, in particular, 
the output as a function of the input knowing the structure and the internal set of parameters 
of the module. It turns out that the dynamics of such RC works as a special kind of cellular 
automata and, due to a celebrated result in [15], they can emulate a Universal Turing Machine! 

For the FC modules, we give a complete description of the dynamics for the cases of one and 
two units in the fundamental case of constant inputs and for all the parameter values of the 
circuits. Using the exponential contracting rate of the system, the behavior for more general 
input signals follows. Here the input level fixes the dynamics of the module in a subset of the 
phase space allowing dynamical regimes that do not correspond to the typical autonomous 
version. 

The relation between the transducer like point of view used in the RC analysis and the 
constrained phase space used for the FC is done by symbolic dynamics. 

It is worth to notice that, as we shall see, modules that we show to perform different dynam- 
ical behaviors may be functionals in different situations [IB], therefore allowing to perform 
different functions either for given values of parameters through different inputs, or vice- 
versa. 

The paper is organized as follows: section 2 recalls the basic properties of the models, section 
3 describes the dynamics of the regulatory cascades (RC), section 4 describes the dynamics 
of the forced circuits (FC) and section 5 shows how, in simple cases, it is possible to treat a 
mixed type module. Section 6 is concerned with final comments and outlooks. 



2. Description of the Model and General Properties 
2.1. Discrete time regulatory networks. 

We consider a special class of models, i. e. discrete-time regulatory networks. These are 
discrete-time dynamical systems on a network (see [El l25j [T4], [6] ) . By a network we mean 
a digraph, with vertices in a given finite set V and with arrows (oriented links) taken from 
another given finite set A C V x V, together with a set of additional characteristics we 
introduce next. 

Vertices account for interacting units carrying a certain activity level (a scalar in [0, 1] asso- 
ciated to each unit v G V), depicting the product of a gene, and the arrows account for the 
interactions between them. 

Although the results presented below hold for more general interactions (as for instance 
multiplicative interactions) we restrict ourselves to additive inhibitory-activating interactions 
for the sake of simplicity (see [25J [HI [6] ) . 
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The model is defined as follows: to each arrow (u, v) G A we associate an interaction threshold 
T uv G [0, 1], a sign a uv G {—1, 1} indicating whether the action of u over v is an activation 
(o~ uv = 1) or is an inhibition (a uv = — 1), and a coupling strength k uv G [0, 1]. The activity 
level of the network at time t G N is specified by the collection of the x* G [0, 1], expressing 
the activity of each unit v G V at time t G N. 

Denoting I(v) := {u G V : (u, f ) G A} the set of vertices acting over a vertex v G V, the 
activation at time t + 1 of the unit w G V is given by 

(1) x< +1 := ax* + (1 - a) Kut ,# K„(x* - T w )). 

The constant a G [0, 1) appearing in the equation plays the role of a degradation rate. In 
absence of interaction, the activity level of a unit v decreases exponentially fast to zero, 
x* = a*x°. In the interaction terms, H represents the Heaviside (step) function, H(x) = 1 
for x > and H(x) = for x ^ 0. Hence each interaction term is a piecewise constant 
function whose value changes whenever one of the coordinates x^ crosses its own threshold 
T 

J- uv • 

Without loss of generality the coupling strengths are normalized by YluelM K uv = 1 for each 
v G V. 

In [6] we described the dynamics of some of this networks in great details. We shall come 
back to these results when comparing with their open versions. 

2.2. Open networks. 

As mentioned above an open network is defined as part of a larger regulatory network together 
with the corresponding incoming and outgoing arrows, respectively from and to the rest of the 
network. Motifs are such open networks that appear more often than the expected frequency 
in some null statistical graph model |17l [3]. In the present paper we study open subnetworks 
without any reference to their possible abundance inside a larger regulatory network and we 
rather call "module" any ('simple') open subnetwork. 

Any collection of N vertices V^od £ V of a network V defines an open network (as long as 
not all the arrows with head in Vmod have tail in V„iod)- 

Formally an open regulatory network or a module is a regulatory network with vertices V mo d £ 
V, the internal units, and three kind of distinguished arrows: 

(1) the incoming arrows, denoted A- ia is the set {(u,v) : u G V \ V mo &-, v G Kood}; 

(2) the inner arrows denoted ^4 mo d is the set {(n, v) : u G V mo( ±, v G Kaod} and, 

(3) the outgoing arrows denoted A ont the set {(n, v) : u G V mo( ± , v G V \ V mo d}- 

The dynamics of the module is defined using the same rule (see Eq. [1]) as for the general case. 
However, as we consider the module as an open system (as if there were no feedback loops 
from V^d to V^d through V cxt := V \ V mo d) we may replace the exact knowledge of x^ for 
u G V cxt by the one of the corresponding symbol: Q uv '■ — -ff(o"ui)(x^ — T uv )). 

Therefore, the state of the system at a given time t G N is determined by a vector x* G 
[0, l]# v m°d and the imposed external activation levels Q l uv for all (u, v) G Ai n . We refer to the 
array of external activations, {9 UV for all (u,v) G Ai n }, as the input code 9- m . 
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In the case of an open system, the evolution rule ([I]) reads: 
(2) x^ +1 :=ax* + (1 - a)A, (x*, 4) 

with 



(3) 



uei(v)nV mod uG/('u)nv r ext 



The state of the internal units determine the activation of the internal arrows A mo( j. 
We put together this data in the internal code 

(4) e l mod := (H(a uv (xi - T uv )) : («, v) E ^ mod ) . 

Finally, the influence of the open subnetwork (y mo &, ^4 mo d) over hs environment is codified 
in the sequences of output activations, which we group in the output code 

(5) dt == (^(cr™(x* u - T uv )) : (u, v) € A Vmod ^ Vcxt ) . 

In Figure [J we give an example of a module together with its external complement. 




Input code 




Input code 



Output code 



FIGURE 1 . Open subnetwork (green) extracted from a larger one (red, green and 
blue). Hammer-like arrows represent inhibitory interactions, standard arrows activa- 
tory ones. 



From Equations ([2]) and ([3]), we see that each possible input code 9 := m uniquely determines 
an affine contraction Fq : [0, l]#^m°d _ > [0, i]#Knod given by 

(6) F„(x) =ax+(l-a)A,(x,0 in ). 

The collection of all these affine contractions defines an iterated function system (IFS) 

(7) T := {Fg : [0, l]# v ^ [0, l]#^°d . q £ jo, 1}# A ™ } , 
from where the dynamics reads: 

(8) x' +1 = F e Mo...oF 9 „(x°), 



6 



R. LIMA, A. MEYRONEINC & E. UGALDE 



with attractor 



(9) 




Fgt-l o 



■■■oF fl o([0,l]) 



Remark 2.1. The form in Equation ([3]) is not the most general for which the results below 
can be proved. Notice however that by a proper choice of the internal and external thresholds 
it already includes the cases of the AND and OR logical outputs [3]. On the other hand, the 
essential ingredient needed for all the proofs is the collection of affine contractions indexed 
by the possible forcing codes, defining an iterated function system (IFS) as in ([7]) and ([8]). It 
is also well known that such IFS are skew product dynamical systems [4]. However, because 
in our case we are interested in properties of very specific models, known results about the 
latter systems are only used to fix the general context of our study. They are, in particular, 
implicitly used when arguing on the genericity of some properties proved below. 

In fact, for discrete-time regulatory networks, there is a way to understand the dynamics of 
an open system in terms of a collection of autonomous mirror systems. For an open system 
(V mo d, A mo d, Ai n ) and any constant forcing code: 9° ■ • • 9 t ■ ■ ■ , with 9 l = 9 Mt € N, let fig be 
the corresponding attractor, named basic attractor and defined by Equation (|9|). Now, for 
each constant forcing code, the system will evolve inside an invariant subset of its phase space 
as an autonomous system, up to an affine change of variables. This invariant subset attracts 
all the trajectories starting outside. The lasts correspond, by the same change of variables, 
to those starting outside [0, 1] in the closed system and for which the behavior is known [6]. 

On the other hand, the attractors for the closed systems, and therefore for each constant 
forcing of the open system, are unions of the so called global orbits, those that can be extended 
backward in time up to — oo inside [0, \\# Vm ° d (see [6]). The same is then true for every basic 
attractor of the open system. The generic case, in measure sense in parameter space, is such 
that each of the basic attractors is uniformly bounded away from the discontinuities (the 
internal thresholds). In this case, any orbit with constant forcing will approach the attractor 
exponentially fast, with rate logo. Therefore, for a general forcing the corresponding orbit 
will wander around the basic attractors and will closely approach one if the forcing input 
stays constant during a sufficiently long duration, and the dynamical behavior of the system 
is clear in this case. Moreover it tells us that a is the main parameter controlling how fast the 
system respond to a new input signal. This suggest that a may be experimentally estimated 
when the corresponding proteins are not actively degraded, as it is the case for most proteins 
in growing bacterial cells [3]. 

For the remainder exceptional cases of the parameter values, the attractor will be arbitrarily 
close to the discontinuities. In this situation an orbit approaching a basic attractor will 
sometimes be close but on the opposite side of the attractor. This will cause an accident in 
the internal code followed by an unforeseeable length of time before the orbit to approach 
again the basic attractor. 

The rigorous description of the dynamics of all the possible internal codes of the attractor in 
this case is still an open mathematical problem. 
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Recall that a path in a network (V,A) is a sequence of vertices (vq,vi, . . . ,V£-\), such that 
(vi, v i+i) € A for ^ i ^ I — 2. If in addition (ve-i, vo) € A, the we said that (vq, v±, . . . , V£—i) 
is a cycle. Circuits are cycles with no repeating vertices. 



3. Regulatory cascades 
3.1. Description of the regulatory cascade. 

Among the open structures encountered in the analysis of biological networks, and that 
could be associated to certain functions, are the feed-forward loops, the dense overlapping 
regulons, and the diamonds (see [3]). These motifs have a two-common feature once regarded 
as digraphs: (a) they are connected and (b) they do not contain cycles. As mentioned in the 
introduction, and as we will show below, open sub-systems on these kind of networks behave 
as finite state translators, from the input to the output code, with a certain delay which in 
principle may depend on the input structure. 

A regulatory cascade (RC) is an open subsystem (V mo d, Anod) defined in a regulatory network, 
such that it is connected and with no cycles. It is clear that such a network cannot be strongly 
connected. The open subnetwork in Figure [1] is an example. 

The simplest RC (see Figure [2]) consists of a single internal vertex v £ V, forced by ver- 
tices U := {ui, . . . ,u n } C V cxt = V \ {v}, and affecting a collection of vertices W := 
{wi, . . . , w m } C V cxt - Its functioning can be though as the translation of an input code 6f n := 
(i7(cr u „(x^ — T uv )) : u 6 U), into an output code #* ut := [H{a vw {^ v — T vw )) : w G W). We 
assume that the vertex v does not interact with itself, so that (vo,vo) £ A. We will refer to 
this open subsystem as as the elementary transducer (ET). 




FIGURE 2. The elementary transducer: an open subsystem consisting of a single 
vertex with several input and output arrows. Lines with hammer like heads represent 
inhibitory interactions, arrows activatory ones. 
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Another elementary regulatory cascade is what we call a regulatory chain RCh. It is defined 
over a linear network — > v\ —>■■■—> v n — > consisting of a path connecting the input vertex 
(t>i) to the output one (v n ), with input and output arrows, (u, v±) and (v n ,w) respectively. 
Since to each interaction mode (activation/inhibition) it corresponds a sign (+1/-1), we can 
therefore associate the sign a : = n&=i °~v h v k+1 x o~ VnW to the RCh — > v\ — > • • • — ► u n — >. As 
we will show below, the functioning of this chain as a transducer essentially depends of its 
sign. 

A remarkable family of RCs are the so called feedforward loops (FFL), consisting of three 
internal vertices connected as indicated in Figure El There are single input and output 
arrows, and two chains (elementary paths) connecting them. Taking this into account, the 
feedforward loops were classified as coherent and incoherent, depending on whether the sign 
of the chains composing them are of the same or opposite signs (see [3] for details). 




• v 



#out T 

FIGURE 3. Incoherent feedforward loop of the type 3 



As we will show below, it is convenient to decompose a given regulatory cascade into elemen- 
tary transducers and regulatory chains. 

3.2. Dynamics of the regulatory cascade. 

We investigate how efficiently and robustly a regulatory cascade can work when transmitting 
information. Two mathematical questions are clearly related to this: what is the time length 
of the external signal needed by the cascade system to tie an output, and what is the possible 
dependence of the output on the internal state of the cascade at the time when the input is 
detected. In particular, our aim in this section is to determine the conditions under which 
each output code of a RC depends only on the sequence of input codes received during a 
finite period of time, whose length depends on the structure and on the parameters of the 
RC. When this phenomenon occurs we say that the cascade resolves the output code with a 
finite delay. To this aim, we will first analyze in details the simplest RC, i. e. the elementary 
transducer. Since every RC can be decomposed into a collection of ET's, the conditions under 
which a general regulatory cascade resolves the code, and is functioning as a code transducer, 
will be deduced from the behavior of its elementary transducers. 
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3.2.1. The elementary transducer. 

For the ET in Figure [21 a temporal sequence of input codes Gi n := (Of n )teN, completely 
determines the evolution of the internal vertex: 

(10) x^ +1 := ax* + (1 - a) E K uv 9 f uv . 

Hence, the temporal sequence of output code 6 out := (0l ut )t&N can be computed from this 
sequences of input codes and the initial condition x°. 

When the input signal has been present for an infinite time, the output code depends only 
on the infinite sequence of input signals. Otherwise, if the input signal started its action at 
time to, then the output code will also depend on x*°, the activity level of the internal vertex 
at time to. 

According to Equation ([10]) . from a sequence of input codes 0j n := (^)* =to , and an initial 
condition x*° , we obtain 

(11) x* = o*-* ^ + (1 - a) a'""" 1 ( E > 

T=t \u£U J 

from which we compute the output code 0* ut := \H [a vw (x* — T vw )) : (v, w) G A out ) • Hence, 
in order to compute the temporal sequence of output codes 6 out := (0* ut )t e N, we have to 
determine the position of x* with respect to the output thresholds {T vw : (v,w) G A)ut} ; 
and so for each time t G N. To solve this problem we take into account the following. 

1. - The internal vertex IFS. Each possible input code 9{ n := (9 UV : (u, v) G ^4i n ) ? uniquely 
determines an affine contraction Fg : [0, 1] — > [0, 1] given by 

(12) F e (x) = ax + (1 - a) ^ k uv 9 uv . 

{u,v)eA- m 

As already mentioned, the collection of all these affine contractions defines the iterated func- 
tion system (IFS) 

(13) :F:={F fl :[0,l]-[0,l]: 9 € {0,1}***} , 
with attractor 

(14) ^:=H |J Vo...of„([0,l]). 

t>0 6)0...6)t-l 6 ({ ,l}#^m)< 

2. - Resolution of the output code. For each input sequence (9{ n ) T ^o an d a U time t ^ 1, 
let Ino ...gt-i '■= Fnt-i o • • • o Fqo ([0, 1]). The evolution of x* is such that for each to G N fixed 

in' in in in 

and all t ^ to we have 

in m "i n "i n 

where 

(15) x,(C t0 • • • C 1 ) == (1 - a) E at " T_1 ( E **C J • 

r=t-to \ueu / 
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the to approximation to the activity level. From Equation (|15|) it readily follows that 

\xt-x v (el t0 ---e^ 1 )\^a t °, 

for all t > to- Hence, for each to fixed, the output code can be resolved for all the input codes 
C*° ■ ■ ■ C 1 satisfying 

mm U v (O t - t °---0 t i ~ 1 )-T uw \ >a^. 

As mentioned above, in this case we say that the code is resolved with a finite delay to- 
The set of input code sequences which can be resolved with a delay to grows exponentially 
with to- Depending on the parameters of the system, the complement of this set could grow 
exponentially as well. 

There are two cases where we can show that the ET resolves the input code in a finite time. 
The first case relies on an internal characteristic of the ET we call internal separability. In 
this case, with probability 1, the thresholds are separated from the attractor. In the second 
case we assume that the input signal satisfies a property we name low input complexity. 

3.- Internal separability. A possible simplification occurs when the attractor is a Can- 
tor set. This is the case if in Equation (| 12f) the contraction rate a is sufficiently small. In this 
cantorian case, for any T ^ fijr there exists a depth tx € N such that T ^ lj0°---0-'"- 1 e({o i}#A in ) T I g o ...q^- 1 -, 
for each r ^ t-r- Hence, if {T vw : (v,w) G A out } (iQf = 0, which happens with probability 1, 
the maximal depth to := maxy tr is such that #* ut = (H (a vw (x v (9 t i ~ t ° ■ ■ ■ 0\~ ))) : (v , w) € A out ) . 
In this case the elementary transducer acts as a cellular automata, transforming sequences 
of input codes to sequences of output ones, with a delay to- Indeed, we can define <& : 

({0, l}#^n) N ^ ({Q 5 l}#^out) N guch that 

(16) mn)r ■= {H (a vw (x^-* • • • C 1 ))) : ™ e W) = Ct*°, 

with x 1 ,(0j~*° • ■ -^j*" 1 ) the to approximation to the activity level defined in Equation (|15p . 
Let us emphasize that in the cantorian case, the condition {T vw : (v,w) £ j4 ut} H ^jf = 
holds with probability 1 with respect to Lebesgue, i. e. a cantorian ET typically operates 
as a cellular automata. If so is the case, we say that the ET satisfies internal separability. 
If on the contrary {T vw : w 6 W} f) Q?r 0, there would be for each t £ N input codes 
(^in ' ' ' ^in" 1 ) that cannot be resolved in finite time. 

Example 3.1. Let us illustrate the functioning of the elementary transducer in the internally 
separable case. Consider an ET with two input arrows (ui,v) and (v>2,v), with coupling 
constants k UiV = 1 — k U2V = 2/3, and a single output arrow (v,w), as shown in FigureHl Let 
T vw = 83/150 and o~ vw = —1. Fix the contraction rate a = 1/5. In this case the associated 
IFS is 

T:= {Fqo(x) =x/5,F i(x) = x/5 + 4/15, F 10 (x) = x/5 + 8/15, F n (x) =x/5 + 4/5}. 

Its attractor Qj? is a Cantor set of box dimension c4>ox(^:f) = log(4) / log(5). In this particular 
case, the third approximant |Lo g i d 2 G( t Q -n#c/\3 lap a} q2 of the IFS's attractor is the disjoint 

in in in ^ ' ^ ' in in in 

union of 64 closed intervals of length 1/125, indexed by codes in ({0, l} 2 ) 3 . None of these 
64 intervals contain the output threshold, therefore we have a depth to = 3 in this case. 
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The translating cellular automata ^({0, 1} 2 ) N — > {0, 1} N is denned in this case by the local 
function : ({0, l} 2 ) 3 -» {0, 1} such that 



ft 1 f) 2 ) 

in "in "in/ 



1 if 0P G {00, 01} or 0° 0^2, G {(10, 00, 00), (10, 00, 01)}, 
otherwise. 



Here, for instance, all input sequences in {11, 01} N produce the same output sequence 000 • • • . 




4-- Low input complexity. The attractor Vljr is a Cantor set if and only if U eg | 1 i#A in -Pe([0, 1]) C 
[0, 1]. When on the contrary Uee{o i}# A m -^e([0, 1]) covers [0, 1], even if the number of input 
codes that can be resolved with a delay to grows exponentially fast with to, the cardinality 
of its complement could also grow exponentially fast. In this case, in order to resolve the 
output code in finite time, we consider a particular class of input sequences. 

The temporal complexity of a sequence of input codes Gj n = (0f n )t^o is defined as 

(17) Ce ln (t) = #{^ n ---^+* : tGN}. 

If the sequence of input codes Oj n is such that for a fixed k G N and all t sufficiently 
large, Ce in (i) ^ t k (in this case we say that the sequences of input codes have polynomial 
complexity), then the output code can be typically resolved with a finite delay. We are 
interested in the polynomial case due to the fact that strongly contractive regulatory networks 
produce sequences of codes with polynomial complexity [14]. We conjecture that in all cases 
the sequences of codes produced by a regulatory network have polynomial complexity. 

We state our result in this case as follows. For sequences of input codes 6i n satisfying 
Ce in (t) t k for a fixed k £ N and all t sufficiently large, then, with probability 1 in the 
output thresholds, there exists a fixed delay time t\ ^ to such that 

C = {H (*W (x w (t (l • • • O)) ■■ weW), 
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with x 1) (^ 1 < ~ <1 • • • ^j" 1 ) as defined by (fT5j) . The ET typically functions as a cellular automata 
when restricted to sequences of input codes with polynomial complexity. In this case we say 
that the system has low input complexity. 

Remark 3.1. The behavior of any totalistic cellular automaton can be obtained from an in- 
ternally separable ET. Since totalistic cellular automata can simulate any Turing machine 
it follows that ET's have universal computing capabilities. 



3.2.2. The Regulatory Chain (RCh). 

As mentioned above, a regulatory chain is an open regulatory network defined over a linear 
digraph — > v\ v n — connecting the input vertex {v\) to the output one (v n ), with 

input and output arrows (u,V\) and (v n ,w) respectively. To this chain we associate the sign 
a := njc=i a v k v k+1 x &v n w The functioning of this chain as a transducer essentially depends 
of this sign. 

We have two possibilities depending on the common contraction rate o G [0, 1). The simplest 
one occurs when a < 1/2, in which case all the vertices v±, . . . , v n considered as elementary 
transducers typically act as a cellular automata. Indeed, if a < 1/2 to each vertex we associate 
the same dyadic IFS 

(18) ^ dyadic := {F e : [0, 1] - [0, 1] : 9 G {0, 1}} , 
with attractor 

(19) dyadic : = fl U ° • • • ° ^''([0, 1]). 

t^o 0°...0*-ie{o,i}* 
This is a Cantor set with box dimension log(2)/log(a~ 1 ). 

In the typical case, when {T VkVk+1 : 1 ^ k ^ n + ljnSldyadic = 0> we can associate, as before, 
to each internal vertex a depth 

(20) t k := max{i ^ 1 : T VkVk+1 G U e o... e t-i e{OA yF e t-i o ••• oify)([0, 1])} . 

In the present case we can recursively define the tk approximation to the activity level, and 
the internal code as follows: 

t-i 

(21) ^(^-C 1 ) := (I -a) £ a^OL 

T=t-tl 

®v k v k+1 '■= H ((Jv k v k+1 {^-v k {0 Vk J\ Vk ■ ■ ■ ^v k ^ 1 v k ) ~ T Vk _ lVk \\ 

t-\ 

*v k+1 (0v k vlX\ ' ' ' 0v k v k+1 ) := (1 ~~ a ) '^2 a * 1 ^ k v k+1 i 

T=t-t k+1 

for 1 ^ k < n. According to this, the (k + 1)— th internal vertex, considered as an ET, works 
as the cellular automata : {0, 1} N — > {0, 1} N , such that 



[ ] k[ h ~ I (1 + <Tv k+1 v h+2 )/2 otherwise. 



Vk+lfk+2 
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It follows that the output code of the whole RCh can be resolved with a delay to := Ylk=i tk> 
by using the composition <E> := <E> n o $ n ^i o • • • o $ x : {0, 1} N — ► {0, 1} N . 

The structure of each one of the cellular automata : {0, 1} N — > {0, 1} N is such that if 
ft" • • • 0j~ is a constant sequence, then 



(23) *(e) f 



(1-<t)/2 if • • • ^f" 1 = 00 • • • 0, 
(1+<t)/2 if^---C 1 = ll---l, 



where a = YYk=i a v k v k+1 x o"« n «; is the sign of the chain as defined above. 

As for the ET, even if the dyadic attractor Od ya di c is not a Cantor set, the output code of the 
regulatory chain can be typically resolved for sequences of low input complexity. In this case 
the action can also be obtained as the composition $ := $> n o$> n _±o- • -o&i : {0, 1} N — > {0, 1} N . 

3.2.3. The general regulatory cascade. 

As mentioned above a regulatory cascade is defined over a connected digraph with no cycles. 
The vertices of this digraph can be hierarchically organized, so that the input code for vertices 
in the k-th. level are the output code of the vertices in the (k— l)-th level. This way, following 
the same idea as in the analysis of the regulatory chain, we can resolve the output code of the 
cascade by using a composition of cellular automata associated to the vertices of the cascade. 

The hierarchy of vertices is the following. On top we have the root vertices, 

(24) V~ root :={v £ V mod : (u,v) £A^u(£ V~ cxt }, 

which are the vertices in V mo d all of whose incoming arrows have tail in the set of external 
vertices. Let Vq = V TOO t, Uo = V~o, and for each k ^ 1 we define 

(25) Vfe := {v G Kiod : (u, v) G A u G C4_i U V cxt }, 
U k := LVlUFfc. 

We can see Uk-i as the union of all vertices up to the (A; — l)-th level. Then the fe-th level, 
Vk, is composed by those vertices in V mo d all of whose incoming arrows have tail in levels 
lower than the k-th or in V eX t- Since the underlying digraph has no cycles or loops, then 
these levels are nonempty. Also, since the digraph V mo d is a finite set, there is a finite number 
of levels, all of them of finite size. Vertices of the last of theses levels are called leaf vertices. 
This last level can also be defined by 

(26) Vi^ := {v G V mod : (v,w) G A^ w G V cxt }. 

The depth d of the regulatory cascade is the number of steps needed, starting from the root, 
to determine the leaf vertices, i. e. Vd = VJeaf- 

Because of this hierarchical structure, and taking into account Equation ([2]), the activity level 
x* , for v G Vfe, is given by 

(27) x^ +1 := ax* + (1 - a) i k w F(ct w (x* - T uv )) + ^ fc TO ^ 

I «e/(t))ni/ fc _ 1 u€i"(v)nv cx t 
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If, on the other hand, we can resolve the internal codes (0*^ : u G I{v) n Uk-i), then the 
previous equation reduces to 



x* +1 := ax* + (1 - a) < 



2 

k weJ(v)nr/ fe _i 



+ 

ue/(i>)nv ox t 



Hence, all the internal vertices can be considered as ETs, and all what we discussed in 
paragraph 13.2.11 applies. Once again we have the alternative between the internally separable 
and the non-separable case. 

1.- Internal separability for RC. Once again, the simplest case occurs when the IFS 
associated to each one of the internal vertices, 



Fg(x) = ax+ (1 



has a Cantor set £l v as attractor. In this simple case, with probability 1 none of the thresholds 
associated to internal arrows will lie inside its corresponding Cantor set, so that to the 
sequence of output codes from vertex v, (0* , : (v, v') £ A) can be resolved with a finite 
delay t v , by using a cellular automata $„ : ({0, 1} #/ ^) N -» ({0, 1}#°(«)) N . Here I(v) and 
0(v) are respectively the input and output set of the vertex v, and the cellular automata & v 
is defined in the same way as in Equation (|16p . The input sequence of the cellular automata 
associated to the internal vertex v € Vk is obtained from the input sequence of codes via the 
action of the cellular automata $ u associated to vertices in levels lower than k. 



Example 3.2. Let us illustrate how operates a regulatory cascade in the internally separable 
case by considering the incoherent feedforward loop of Figure El This RC contains a positive 
regulatory chain — > u — > v — >, and an negative arrow from u to w. Let us fix a = 1/5, output 
threshold T WVout = 1/2, output sign a WVout = 1, and coupling constants k uw = 1 — k vw = 2/3. 
In this way, the functioning of the vertex w considered as an ET could be deduced from that 
of Example 13.11 

Let T uv = 22/25, T uw = 22/125, and T vw = 1/2. The IFS associated to these vertices is 

T u = Tv = {F Q x = x/5, F lX = x/5 + 4/5}, 

which has a dyadic Cantor attractor Q v with box dimension log(2)/log(5). The delays de- 
termined from the position of the thresholds with respect to the attractor are t u = 2 and 
t v = I- 

The internal code {9 l uv , l uw ) depends only on the input sequence Oj n and it is given by 

t / o if et- 2 el-],G {00,01,10}, 

uv 



i ifC 2 C _1 = ii- 

it— 3 at— 2 nt— 1 



1 if C C C e i 000 ' ° 01 ' 01 °}' 



uw I otherwise. 
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FIGURE 5. Behavior of the incoherent feedforward loop. On top, the input signal 6^. 
Below, the state of the output vertex x^ responsible of the output signal. The output 
signal changes each time crosses the output threshold indicated by a horizontal 
line in the figure. 



For the RCh — > u — > v — >, the composition of & v o <3? u allows to determine the internal code 
Q\ w directly from the input code Gj n . We obtain the following: 



'vw 



if e\-*o\- 2 e{oo,oi, io}, 

1 otherwise. 



Finally, considering the vertex w as an ET with two inputs we obtain an output code 

it-int-3 



lif er*^ € {oo,oi,io}, 
o if efflfc 3 = ii- 



This feedforward loop operates as a cellular automata. The output code can be resolved with 
a delay t w = 4. 



2.- The general case. The output code cannot be determined solely from the input code 
when at least one of the attractors £l v is not a Cantor set, or when one of the thresholds T vv i 
is contained in its corresponding attr actor Vl v . In those cases, for each depth r € N there are 
sequences of input codes such that the output code cannot be resolved in finite time. The 
comment made in Paragraph 13.2.11 concerning sequences of input codes with low temporal 
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complexity applies once again. In that case we can typically resolve the output code in finite 
time. 

Remark 3.2 (Comments on decrease of complexity). Each time we can resolve the output 
code in finite time, whether we are in the internally separable case or because the low input 
complexity holds, the RC operates as a cellular automata <E> : ({0, l}# A in) N — > ({0, l}#^°ut) N 
such that 

The symbolic complexity of the output code is the counting function Ce out : N — > N such 
that 

(28) C 0out (i) := # {Ct • • • C* : T € N ) • 

i. e., for each t £ N, Ce out (*) counts all different sequences of output codes of length t. In 
the present case, since the output code can be resolved with a delay to, we clearly have 

(29) Ce out (t + t )^Ce in (t) + (2* w ) t0 . 

In the case of an exponentially increasing complexity the previous inequality ensures the 
non-increase of the entropy 



4. Forced Circuits 



Forced circuits (FC) are modules in a regulatory network such that the underlying graph 
(V mod ,A mod ) is a circuit. 

Examples of forced circuits are the (open) Negative Auto Regulator (NAR), the Positive Auto 
Regulator (PAR) and the network motifs with double-positive (or double-negative)-feedback 
loop [3]. 

In this section we focus on special cases of FC in order to show the strategy of analysis as well 
as the reach possibilities of dynamical behaviors displayed by such modules. Not surprisingly, 
each one may show different behaviors according to the inputs. In particular, in most cases 
the dynamical behavior of the module (open system) is different from that of the same system 
when isolated. 

4.1. The self regulations. 

Self-regulation occurs when a transcription factor acts as an inhibitor (self-inhibitor or NAR) 
or an enhancer (self-activation or PAR) of the transcription of its own gene. Self-regulation is 
a very common situation, for instance it is involved in over 40% of known E. coli transcription 
factors |22j . 

For a self-regulator subjected to a single input arrow as indicated in Figure El thanks to the 
normalization condition k vv + k uv = 1, the Equations ([I]) and ([2]) read: 

(30) x< +1 = ax* + (1 - a) [H (<r w (x* - T m )) + K uv (0^ - H (<j to (x* - T w )))] , 
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From (JTJ) , we see that for re. 



one recovers the isolated self-regulation. 



£3 



9 



FIGURE 6. The open self-regulation, left: a vv = —1 for a self-inhibition; right: 
<j vv = +1 for a self-activation. 

Otherwise, for each given a vv there are two extra parameters, the input intensity k uv and 
the input signal sequence Bi n := (^)tgN (also called "exogene variable"). As usual we write 
H (a vv (jc v — T vv )) = 9 VV and then we merge internal and external codes in a unique symbol: 
9 t = 9* v 9^ v (internal on the left, input on the right). 

As a first approximation, an easy way to visualize the different dynamics when changing the 
parameters and/or the input signal sequence is to show a diagram with the possible transitions 
among the 9s at each time step (see Figure [9] (B)), known as dynamical graphs \23\ \9\ IT]. In 
particular, it allows to locate forbidden paths in phase space and therefore localize the most 
robust bifurcations. Because residence times in each loop are not specified in this diagram, 
dynamical graphs carry only part of the dynamical information. 

We shall describe in the following the case of the self-inhibition in detail and then, for con- 
ciseness, we only sketch the case of the self-activation. 

4.1.1. The open self inhibition. 

The dynamics of an isolated self-inhibitor, NAR, consists only of oscillations [6]. In fact, 
whenever < T vv < 1, this system is conjugated to a rotation on a circle with a rotation 
number u(a,T) depending on the parameters (5j[6]. 

Figure [7J illustrates the three possible cases of the dynamics for the open self-inhibitor, 
shown Figure EJ when varying the input intensity. Accordingly, the parameter subspace 
€ (0,1), a G [0,1) and k uv G [0,1]} can be divided into three input 
intensity regions, corresponding to different dynamical behaviors. Region I: if k uv <C T vv , 
Region II: if T m , < k uv < 1 — T vv and Region III: if k uv > T vv . 

Remark 4.1. Up to a change of (9 VV , 9 UV ) in (1 — 9 VV , 1 — 9 UV ) in the dynamical graphs, it is 
enough to consider the cases where T vv < 1 — T vv . 



Figure [8] displays the response of the self-inhibitor circuit to an input in Region I and Region 
III. The input 9 l u is set to for ^ t < 20, then equal to 1 for 20 ^ t ^ 40 and then again to 
for 40 < t < 50. Therefore, the dynamics is governed by the branches /oo and /io during 
the first and the last period and by /oi and /n during the intermediate time. It is clear from 
Figure [7J that the response is very different in case (A) of small input intensity and (C) of 
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(A) (B) (C) 

FIGURE 7. The open self-inhibition: graphs of the IFS F 6uv (x v ) = ax t + (1 - 
a) [H (T vv — x.„) + k uv (6 UV — H (T vv — x„))] For, (A) Region I, small input intensity 
k uv , (B) Region II, intermediate input intensity and (C) Region III, high input inten- 
sity. Possible fixed points are given by the intersection of the graph with the diagonal. 
The notation fij,i,j = 0,1 stands for the branch of Fi when H (T vv — x„) = i and 



high input intensity. In the first case there are no available fixed points and the system runs 
in pure oscillations and only the amplitude and the frequency distinguish the lower from the 
upper level of the input. In the last case the system contracts to the lower fixed point (see 
/io) for lower level activation and to the upper fixed point (see /oi) for high level input. It 
is clear in this example that the input intensity may change the dynamical behavior of a 
module from oscillatory to bistable. 
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FIGURE 8. The open self-inhibition: The typical response x* to the input x^: 
(A) in Region I, small input intensity k uv , (C) in Region III, high input intensity. 
The horizontal green line indicates the output threshold. 
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In Figure [U these three regions of parameters are subdivided in smaller subregions corre- 
sponding to different dynamical graphs and therefore to possibly different dynamical regimes 
that we describe in the following. 
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FIGURE 9. For the open self- inhibition, (A) Parameter space. It is made of 13 sub- 
regions. Different colors, depicted for a given T vv , show how the k uv parameter can 
affect the dynamics. (B) Corresponding dynamical graphs of possible transitions. The 
code denotes 9 VV 8 UV . Colours are in correspondence with the plane in (A). Region 
I: includes 1, 2, 6, 7, 11, region II: 3, 4, 8, 9, 12 and region III: 5, 10, 13. 



1 .- Region I. For k uv < T vv , corresponding to Figure [7] (A) , the dynamics is the same as for 
the isolated self-inhibitor, oscillates whatever the input sequence can be, but with 
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an amplitude and a frequency that depend on the input (^)teN- This is because in this case 
the corresponding maps have no fixed points. 

If a is small enough (subregion 1), the system will change code at each time and the next 
code is completely defined by the external input at that time, Figured] (B). In particular, for 
a constant input sequence the coding of the oscillating orbits (x* does not depend on the 
initial condition x°. 

For higher a values (subregions 2, 6, 7 and 11), the input sequence (9^ v )t£^ can non-trivially 
affect the dynamics of x v . This corresponds to the occurrence of the loops in the dynamical 
graph, allowing the orbit to lie inside a given atom of the symbolic partition more than one 
time step depending of the input sequence. The reason is clear: as a increases the image of 
the branches (see Figure [7] (A)) may intersect the two sides of the discontinuity. 



Remark 4.2 (The general case: local fixed points and absorbing intervals). One can consider 
the more general case of any number of inputs. Again the existence of local fixed points for 
the IFS Fg uv depends on the internal and external parameter values. Let P be the number 
of local fixed points for Fg uv . Let S r (v) be any subset (possibly empty) of I(v) such that 
YlueSJv) Kuv > Tvv and let R be the number of such subsets. Let Si(v ) be any subset (possibly 
empty) of I(v) U {v} containing {v} such that YlueSiM k to ^ T vv and L be the number of 
such subsets. It is not difficult to see the relation: P = L + R. Consequently ^ P ^ 2^ zI ^ . 

Similarly, let A be the number of local absorbing intervals C (0, 1) and let S c (v) be any 
subset (possibly empty) of I{v) and S c {v) = S c (v) U {v} such that X^eS c («) Kuv ^ ^ vv < ~ 
Sug5 c (d) Kuv > an d & be the number of such subsets. Then A = C. Consequently ^ A ^ 

Finally one can check that P + A = L + R + C = 2 #I ^ . The two extrem cases are 
P = A = 2# 7 ^) if and only if 1 - k vv ^ T vv < k vv , and P = 2* J ^ A = if and only 

if Kvv *S Tvv < rnhkte/(u) {^u} < - )r Kvv ~\~ max ng/(i;){ K ™} ^ Tvv ^1 Kvv 

Externally induced switches from (projected) fixed point converging regimes to (projected) 
periodic orbit converging regimes, and vice versa, give fairly simple dynamics. It results 
in the concatenation of pieces of (transient) orbits from either regimes. Switches between 
periodic orbit converging regimes may however be quite complicated as we see it next. 

2.- Region III. For the opposite case, i. e. for k uv > 1 — T VV corresponding to Figure [7] (C), 
we notice the occurrence of two fixed points. The first, denoted xo, corresponds to the branch 
/io and therefore an orbit will come close to it by the repeated injection of the corresponding 
input 6 UV = 0. The second one, xi, corresponds to the branch /oi and an orbit will came 
close to it by the repeated injection of the corresponding input 6 UV = 1. Therefore, in this 
case, the internal code 9 VV is a delayed slave of the input code 9 UV provided the forcing is 
permanent enough. 

In this sense the NAR behaves in this region of parameters as the RC described in the previous 
section. 
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In particular, if (9u V )t&t is constant (i-independent), then there exists a t° ^ such that 
d l vv = 1 — Q l uv Vi ^ t , so that the system can be easily driven in one of the two different 
states as if it was bistable. 

Notice that in the dynamical graphs of Figure [9] (B), for the corresponding subregions 5, 
10 and 11, whatever the input sequence is, the loops that do not correspond to codes 10 
or 01 cannot consecutively be taken an infinite number of times and the (finite) number of 
residence steps will depend on the input sequence after some delay. 

If the input sequence (6*^)t G N is not constant, then the internal code (#£,y)teN will depend on 
the initial condition xJJ. In this case, depending of the parameters, the proper mathematical 
study of the internal code sequences that correspond to a given input sequence is still an 
open problem. 

The converse statement can however be formulated and is also of practical interest. Namely, 
being given an observed internal sequence, what are the possible inputs and internal param- 
eter values of the self-inhibition that realize that observed sequence? 

We now present an example showing that the use of admissibility conditions as in [6], im- 
plemented in a numerical algorithm, allows to produce such sequences for each particular 
occurrence. 



Example 4.1. Let us illustrate the identification procedure in the case of a periodic orbit of 
an open self-inhibition with the observed internal sequence (01001) 00 of period 5. 

First, we indentify within Figure[9](B) the candidate dynamical graphs to realize the observed 
sequence, and then the families of input sequence provided by those graphs. In our example 
the candidate input sequences [Q t uv )t<^'L are the families (with uj G {0, 1}): 

• (ujuj1ujuj) c 

and 



° for 
° for 



(0wl0w) c 
(01101) 00 for _ 
(wwwww) 00 for [6], \j} and |U| , 
(OloujOuj) 00 for [j], and [l2~| , 
(01^01) c 



for 10 and 13 



Notice that no transitions in the graph [T] can produce the observed internal code. 

Then we choose one input sequence among the more robusts. Those are the one appearing 
in the greatest number of dynamical graphs. This criterion makes more likely the set of 
parameter values that realize the observed internal code for that candidate input code to be 
broad in the parameter space. This way we select the candidate sequence (01101)°° common 
to all transition graphs except [T"| , and for comparison we also consider (00111)°° that is only 
possible for |T|, [f[], \j} and [U| . 

Thirdly, we explicit the admissibility condition with the internal and the candidate input 
codes, and solve it. The admissibility condition writes: 

(31) 



sup x ^ T uu < inf x" 

*eZ:(6>* lu ,e*j6{l}x{0,l} teZ:(0* u ,0*„)e{O}x{O,l} 
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with for a fixed (a, k uv ) £ [0, 1) x [0, 1], 
(32) 



{x* : t € 1} = { 1 



1 



1 



1 



j ^ ^ ^(n+fc)mod5 



+ K ra ^a( p+fc ) mod5 



(n,p)e[0,4] : 



In the last expression I s = {k € [0, 4J | a~ 1 {9 uu ) = and a- k - l {6 uu ) = if s = and o-* -1 ^ 
1 if s = 1}, and E = {k £ [0,4] | ct -1 ^) = and o-^- 1 ^) = 0}, a being the left shift 
map in {0, 1} Z . We have explicitely in our case 1$ = {1,2,4}, I\ = {0,3}, and E = {2,4}. 

The results are shown Figure PTOl with the plot of all the (a, T vv ) values admissibles for a fixed 
value of k uv . As expected from the analysis of the transition graphs, the domain for (00111) 00 
appears smaller than the domain for (01101)°° when varying the k uv parameter. 





0.2 0.4 0.6 0.8 1 0.2 0,4 0,6 0,8 1 



a 

(A) (C) 

FIGURE 10. (a,T) subdomains of admissibility for the periodic input (01101) 00 
with K uv = 0.2 (A) and k uv = 0.4 (B), and for the periodic input (00111)°° with 
k-uv = 0.2 (C) and K uv = 0.4 (D). For (A) and (C) the admissible values are within 
the colored areas while for (B) and (D) the admissible values are in between the 
colored boundaries. The area and boundary colors arc in corrcspondancc. 



One can see the overlap of the domains for a given value of K uv (Compare (A) and (B) for 
instance), meaning that the same observed internal sequence can be observed for the same 
parameter values for different input sequences. It illustrates the robustness to that parameter. 

Also, depending on the input sequence there may have a dependence to the initial condition. 
This is illustrated on Figure [10] where each color corresponds to an attracting periodic orbit 
with its bassin of attraction. Furthermore, although all the attracting periodic orbits may 
have the same period, they differ in particular by their amplitude. Another way to present 
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it is that in pratice different realizations may look similar (frequency) but differ (in the state 
at a given time, in amplitude) depending on the initial condition. 

Finally, notice the qualitative relation between the domains of Figure [10] and those of Figured] 
(A). The sequence (01001) 00 is indeed for the input (01101) 00 out of the tetrahedron [T] as 
indicated by the derivative of the uppen bound of the red domain a a = Figure [10] (A). 



Also The sequence (01001) 00 is indeed for the input (00111)°° out of the tetrahedron 8 . 



Remark 4.3. System identification is closely related to the control issue, as if we can find out 
an input and parameter values that realize a given behavior of the module, then that input, 
with the appropriate parameter values, can also be used to control the module to a desired 
behavior. Furthermore the robustness to small aditional perturbations has been illustrated 
on Figure [10] 



3.- Region II. We first treat the case 1 — k uv < k uv . 

If T vv ^ k uv < 1 — T VV , corresponding to Figure [7] (B), the dynamics shares the characteristics 
of both previous systems. This is because in this parameter region only the branch /01 
(matching the input 9 UV = 1) induces a fixed point, while on the branch /oo the system 
will oscillate in any case. That is to say: if 9 l uv = 1 Vi 6 N, then any orbit converges 
asymptotically to Xi, if 9 l uv = Vi £ N, then the dynamics is as for the pure self-inhibition, 
up to a change of variable that depends on k uv (the same as for region I). 

The region II contains three subregions denoted 3, 4, 8, 9 and 12 in in Figure [9] Again, 
depending on the subregion, some transitions, in presence of a suitable input, are possible or 
not. As in the previous case, and for the same reason, only the loop corresponding to 01 in 
the dynamical graph may be taken an infinite number of consecutive time steps. Again, the 
(finite) number of residence steps in the remainder codes will depend on the input sequence 
after some delay. 

In case that k uv < 1 — k uv , due to the symmetry k uv < — > 1 — k uv we shall have in Region II a 
steady low level state xo < T vv corresponding to the input 9 l uv = and high level oscillation 
corresponding to the input 9 l uv = 1. 

Finally, we see that depending on the parameter values of the circuit and on the different 
inputs the circuit may either oscillations, exhibit a low level or a high level steady state, or 
even be bistable. 

This analysis shows that the range of external signals can have a significant influence on 
the expression dynamics of a self-inhibited gene. Depending on the input signal and/or the 
self-regulation parameters, very different dynamical regimes can exist and satisfy different 
functional 'demands' of the whole regulatory network through the subsequent interactions 
the gene is involved in. 
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Xq J-vv Xi v xi 



(A) (B) (C) 

FIGURE 11. The open self-activation. Graphs of the IFS Fg uv (x„) = ax t + (1 — 
a) [H (x„ — T vt) ) + k uv (6 UV — H (x„ — T^„))] . For (A) Region I, small input intensity. 
(B) Region II, intermediate input intensity and (C) Region III, high input intensity. 
Possible fixed points are given by the intersection of the graph with the diagonal. 
The notation = 0,1 stands for the branch of F$ nv when H (x„ — T vv ) = i and 

@UV 3' 



4.1.2. The open self activation. 

We now proceed with the positive self-regulation (PAR). An isolated PAR is a bistable system 
apart from the (non functional) case where T vv £ [0, 1] [6]. 

Figure [TT] illustrates the three possible cases of the dynamics for the corresponding open 
system, shown Figure [6l when varying the input intensity. Accordingly, as in the case of 
the self-inhibition, it follows that the parameter subspace {(T vv , a, k uv ) : T vv € (0,1), a 6 
[0,1) and k uv S [0,1]} can be divided into three input intensity regions corresponding to 
different dynamical characteristics. The corresponding parameter domains are as for the 
NAR: region I: if k uv < T vv , region II: if T vv < k uv < 1 — T vv and region III: if 

Depending of the input intensity, some or all of the possible fixed points are present: and 
xo corresponding to the input 9 UV = (0 is present for all values of the input intensity) and 
(xi and 1 corresponding the input, 6 UV = 1 (see Figure [TT]) . Note that 1 is present for all 
values of the input intensity. 

Remark 4.4. We will consider only the cases were T vv < 1 — T vv for the same reason than 
for the open self-inhibition. 



1.- Region I. For k uv < T vv , corresponding to Figure [TT] (A), there exists two absorbing 
regions, (0, xi) and (xo, 1), for the orbits (x* ) tg N on each side of the threshold T vv whose size 
is fixed by the input intensity k uv . (Notice that if k uv = then the dynamics is that of the 
isolated self-activation.) 

Therefore, for a constant input sequence the dynamics is bistable and an orbit x* converges 
to either or xq, depending on the initial condition, if 6 l uv = Vt, or to either xi or 1 if 
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If 9 l uv is not constant, then the dynamics of x* is driven by the input sequence, with possible 
oscillations within one of the two absorbing intervals (0,xi) and (xo, 1). 

2. - Region III. For k uv 1 — T vv , corresponding to Figure [11] (C), whatever the value of 
a, if the input sequence is constant then the dynamics has one fixed point, at if 9 UV = 
and at 1 if 9 UV = 1. 

In this input intensity region the open self-inhibition is a slave to the input sequence. It 
exhibits, however, a faster convergence to either fixed point if the input intensity k uv is close 
to T vv , [3]. 

As for the self-inhibition in the same range of parameters, if the input sequence (9^)^^ 
is not constant, then the internal code (O^t^n will depend on the initial condition x° and 
again, using the admissibility conditions as in [6], it is possible to produce such sequences in 
each particular (see Example 14. ip . 

3. - Region II. We start with the case 1 — k uv < k uv . 

If T vv ^ k uv < 1 — T vv (see Figure [TT1 (B)), depending on the input sequence (6'*„)t g N, the 
dynamics can be either monostable, with attracting fixed point 1 if 9 l uv = 1 Vi, bistable with 
two possible attracting fixed points and xo depending on initial conditions if 9 l uv = Vi, or 
showing more complicated (oscillating) dynamics if 9 l vv is not constant. 

As for the open self- inhibition circuit if k uv < 1 — k uv , due to the symmetry k uv < — ► 1 — k uv the 
dynamics in Region II can be either monostable, with attracting fixed point if 9 l uv = Vt, 
or bistable with two possible attracting fixed points 1 and xi depending on initial conditions 
if0L = lVi 

Finally we see that, depending on the parameter values of the circuit and on the different 
inputs, the circuit may either oscillate, or be bistable for both inputs, or to show a low level 
or a high level steady state and bistability. 

This classification of the different dynamical regimes may be refined exactly in the same way 
as for the open self-inhibitor. It is clear from the analysis made for the NAR that crossing in 
the parameter space the same limits as for the self-inhibitor will cause robust bifurcations, 
leading to changes in the corresponding dynamical graph. 

Finally, as for the self-inhibition, this analysis shows that the range of external signals can 
have a significant influence on the expression dynamics of a self-activated gene. Depending 
on the input signal and/or the self-regulation parameters very different dynamical regimes 
can exist. 



4.2. The open negative 2 circuit. 

The dynamics of the negative 2-circuit is studied in detail in [6] . It has a very reach dynamics 
including periodic and quasi-periodic attractors that may coexist for certain parameter values 
and there are no fixed points. These attractors organize the dynamics in the phase space. 
We shall describe the main features of the dynamics of the open negative 2-circuit shown in 
Figure [T2"| as compared with the autonomous system. 
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FIGURE 12. The open negative 2-circuit. 

According to the last comment in section 12.21 it is important to understand first the cases 
where the forcing sequences are constant, 6\ n = or 9\ n = 1 for all t. For 9\ n = and 1, 
consider the two maps Fg as in ([6]) and denote x# their fixed points different from and 1. 

There are two completely different cases for that system. 

First, the oscillatory induced regime that works as for the autonomous negative circuit: if 
T m < xo and for as input, any trajectory of the x„ component of Fq will end up after 
a finite number of time steps inside the invariant interval [0, xo) and from there, up to an 
affine change of variables, the dynamics of Fq behaves as the negative autonomous 2-circuit 
studied in [6]: i. e. it oscillates. The symmetric situation occurs for xi < T vv and 1 as input. 
After a finite number of steps the dynamics of the x„ component of Fx oscillates inside the 
interval (xi,l]. In short, the two variables x n and x,j oscillate out of phase (by vr/2). In 
case of forcing they oscillate at a low level and in case of 1 as input at a high level. This 
situation occurs after some transient time if the initial condition of x„ happen to be outside 
the corresponding invariant interval. 

Second, the fixed induced regime which is different from the autonous case: if T vv > xo and 
for as input, any trajectory will be attracted by the unique fixed point of Fq. This fixed 
point is (l,xo). The reason is that, after a finite time we get x* < T vv and from then the 
signal 9 l vu sent to node u is always 1 (since this interaction is an inhibition). As consequences, 
x u — ► 1 and the signal Q l uv sent to node v is then always 1 (since now this interaction is an 
activation). Therefore, also x„ — > xo. Finally, in this case, for a as input, the circuit ends 
up in x u = 1, x v = Xo- The case xi > T vv and 1 as input is solved in the same way, and the 
circuit ends up in x u = 0, x„ = xi. 

Notice that in the fixed induced case, contrarily to the autonomous circuit, it is not the initial 
conditions that determines the final destination of the system but the external input. 

Now, by a convenient choice of the external versus internal intensity it is possible to set 
xo < xi or xo > xi. Therefore there are four possible open negative 2-circuits: (1) the 
bi-oscillating circuit that oscillates at small amplitude for a low level input and at a high 
amplitude for a high level input; (2) the oscillating-fixed circuit that oscillates at low am- 
plitude for low level input and converges to x M = 0, x„ = xi for high level input; (3) the 
fixed-oscillating circuit converging to x u = 1, x„ = xq for low level input and oscillating for 
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high level input and (4) the bistable circuit that converges to x u = 1, x t , = xo or to x n = 0, 
x„ = xi according to a low or high input level. 

Finally, for a more general (time variable) input this open circuit may work in different 
manners and, by appropriately tuning the forcing it is possible to switch the circuit from one 
regime to another. 

By the same arguments, the open negative 2-circuit with reversed sign of interactions works 
in the same manner. 

4.3. The open positive 2 circuit. 

The same type of argument shows that for a positive 2-circuit and a constant external forcing, 
in one case it is bistable as is the autonomous system [6], but becomes monostable in the 
opposite case. Therefore, it is also possible to build bistable-bistable, bistable-monostablc, 
monostable-bistable and monostable-monostable circuits working in the pointed regimes for 
low-high input levels. 

Again it is clear in this case that an open positive circuit may operate in a regime different 
from the autonomous counterpart. In particular, an oscillatory input may drive the positive 
circuit in oscillations. 

Notice that in [1] the same discrete-time piecewise-afhne model have been numerically studied 
to account for the dynamics of a p53-Mdm2 genetic regulatory network made of an open 
positive and an open negative 2-circuit in interaction, with one inward regulation acting on 
each of the two circuits. In our notation, node u stands for p53 and node v for Mdm2. 

The analysis above extend to the case of open n-circuits, in the same way as for the corre- 
sponding autonomous circuits [6]. 

5. Example of mixed type 

In the more general case when a module owns one or more circuits inside a regulatory cascade 
it is in principle possible to combine our results and describe all its possible dynamical 
behaviors. For large networks this method will soon end up in the study of the network as a 
all and therefore it is of any help. However for simple networks the splitting of the all network 
in elementary RC and FC modules may help the study of their dynamics. As an example let 
us consider the incoherent type 1 feeforward loop (Il-FFL) [3], depicted in Figure LT3l We 
first determine the possible dynamical regimes for the open self-inhibitor present at the level 
of GalS. Since there are two external inputs for this circuit (CRP and galactose) its dynamics 
displays four branches, corresponding to the possible combinations of the two inputs. 



Therefore we can organize this information in a 2 x 2 table where say, arrows correspond to 
9cRP,Gais = or 1 and columns to 6 ga iactose = or 1 (notice that, because galactose repress 
GalS, according to our notation, 6 galactose = stands for high galactose level and 9 galactose = 1 
for low level). Now by looking for intensities of the interactions and thresholds, each entry of 
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GalS 



galETK 



FIGURE 13. Incoherent feedforward loop of the type 1 from [3J pag. 454]. 



this table can be filled on the basis of the 3 possibilities described in 14.1.11 steady low level, 
steady high level or oscillation. In each case this information immediately fixes & GalS galETK 
where T2 is the delay due to the interaction GalS, galETK as described in 13.2.31 Notice also 
that CRPGals ls J us t e Q ua l to the input signal O c ^mp after the delay t\ due to the node 
CRP. 

Of course the delays t\ and T2 are only important during a transient time after one of the 
external signals eventually changed (as it is the case for a pulse input). If not, time translation 
invariance of the inputs simplify the analysis. 

It is then sufficient to incorporate this information in the description of the dynamics of the 
RC made in 13,2.31 As a consequence this module can end up with, either an oscillatory 
regime for the final gene operon galETK corresponding to an input corresponding to an 
"oscillation" entry of the table, or in a "steady state" (low or high level) for an input with 
this entry. 

We emphasize that this module may operate in very different regimes either by a change of 
parameters (due for instance to a mutation) or by some modulation of the input signals, or 
both. 



6. Final comments 



In this work we made an attempt to understand how a small regulatory network is operating 
under external stimulus. This stimulus can be, either the output of a larger network where 
the module is inserted or simply an external signal acting as a trigger mechanism for the 
action of the module. 

For Regulatory Cascades (RC) we have shown that there are two (non exclusive) conditions 
under which the RC acts as a finite time delayed transducer (or translator) exactly as a 
cellular automata. It takes a finite input "message" in another finite output "message". 
This two conditions are, either an internal property of the module (readable in its parameter 
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values) or an external characteristic of the input signal. In a sense we can say that under 
one (or both) of this conditions the module works in a safety operating regime since it can 
transmit an unambiguous command with only some delay. 

In the opposite case, the situation is different since then, in principle, the module needs a 
variable, eventually infinite, time to "understand" the input message and then to be able 
to produce an unambiguous output. Moreover this time delay may depend of the particular 
input signal (see final part of Remark 12. ip . No doubt that in this case the module is not 
of a great help for the system, even if this behavior is fascinating from the point of view of 
nonlinear dynamics. But maybe Nature has found already somewhere a situation where this 
fuzzy-like operating regime has some selective advantage! 

From a mathematical point of view our strategy was first to prove these properties on the 
simplest possible module (the elementary transducer ET) and then extend them step by step 
up to the general case (RC). 

For the Forced Circuits (FC) we use a slightly different strategy. In each case we work 
out the subsets of the phase space involved in the dynamics under each external input. 
Here the case of stationary inputs enable a simple classification of the different "extreme" 
dynamical responses of the module. Then thanks to the contracting (or diffusive) properties 
of the dynamics (a consequence of degradation) we may predict the dynamics under more 
general inputs. This also opens the possibility to build small modules with different desired 
functionalities. From a mathematical point of view our strategy was simply to restrict the 
dynamics in each FC to the part of the phase space that became invariant and attractive 
under each specific input. 

These small open circuits show how their functioning may depend on the parameters as 
well as on the different type of inputs. In this context, the notion of modularity for the 
accomplishment of a function cannot be reduced to a simple decomposition of the network 
into subsets of nodes and interactions. 

An example where such ambiguity may be interesting for a biological function is the case 
where for a while a cell population needs two types of differentiated cells (corresponding to 
low and medium expression level of genes for instance, see Figured!] (B), up to a moment 
when, under an appropriate signal (generated by a stress for instance), the two states switch 
to the same effective differentiated state (high expression level of genes for instance), see 
FigureE] (C). 

Now we may wonder why we need one point of view (transducer) when dealing with RC and 
another (constrained phase space) to study FC. The two points of view are in fact equivalent. 
The bridge between them is symbolic dynamics. The starting point of symbolic dynamics is to 
encode the atoms of a suitable partition of the phase space by a set of symbols and to transfer 
the description of the dynamics in the phase space in terms of rules on the corresponding 
symbolic representation. In our case the partition of phase space is naturally determined by 
the set of thresholds and the corresponding symbols by the values of the Heaviside function 
H (<t(x — T)). With this dictionary in mind it is clear that the transducer in the RC case uses 
the phase space made available by the input signal to encode the message and, vice-versa, 
in the FC each node send to the others a code defined by its position in the available phase 
space. These are just two faces of the same coin. 
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Our conclusion is that the existence of open designed dynamical modules leads to an extremely 
reach set of different dynamical behaviors making such units capable of carrying through 
various performances in response to different external stimuli. But, on the other hand we 
have seen that in many cases the same functionality may be performed by different modules. 
The possible criteria for the selection of a particular module among other displaying the same 
function open to very interesting questions [18j . 

The fine tune of the dynamics of a module by a complex input is a current subject of our 
efforts. As for further research in a connected direction, we mention the interesting problem 
of reverse engineering, control and bifurcations for such modules that we have just touched 
in Example 14.11 . 

We are still far from a complete understanding of the way how module's architecture, param- 
eters and inputs are related to functionalities in all cases, and vice-versa. Nevertheless we 
are convinced that our results give an insight in the fascinating interplay between modularity 
and dynamics. 
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